function [yobsCounter,statesCounter]=...
    kfilterRegCounterfactual(parvec,funcmod,Y,...
    solveopt,addsol,etaOrig,aZero)
% Compared with kfilter_full.m, delivered identical results, Jan 9 2012. 
% ====================================================================
% ====================================================================
%% 1. Model solution
% Matrices of second sample stored in structure second 
[G,R,C,eu,SDX,Z]=feval(funcmod,parvec,solveopt,addsol);
% A. Non-existence in either sample 
if isequal(eu,[1;1])==0 
    return
end
%% 3.1. Dimensions and storage 
[NObs,NZ]=size(Y); 
NS=size(G,1); 
NX=size(SDX,1); 

[yobsCounter,statesCounter]=feval(@kfilterRegSplitSimulation,aZero,etaOrig');
%% Subroutine kfilterRegSplitSimulation allows to simulate the model  
% Inputs 
% sInitial: [NS 1] Initial State 
% innovMat: [Nx Nobx] matrix of innovations 
% By being a sub-routine it has access to all variables defined above 
% Be-careful not to use an index (ii,jj) used above or to repeat variable
% names
    function [ySim,sSim]=kfilterRegSplitSimulation(sInitial,innovMat) 
        if ~isequal(size(innovMat),[NX NObs]) 
            error('Input innovMat must be [NX NObs]') 
        end 
        sSim=zeros(NS,NObs);
        ySim=zeros(NZ,NObs);
        sSim(:,1)=G*sInitial+R*innovMat(:,1);
        ySim(:,1)=Z*(sSim(:,1)+C); 
        for kk=2:NObs;
            sSim(:,kk)=G*sSim(:,kk-1)+R*innovMat(:,kk);
            ySim(:,kk)=Z*(sSim(:,kk)+C);  
        end
        ySim=ySim'; 
        sSim=sSim'; 
    end
%% End of File
end 